SPECTRAL METHODS FOR PARTIAL DIFFERENTIAL 

EQUATIONS IN IRREGULAR DOMAINS: 
THE SPECTRAL SMOOTHED BOUNDARY METHOD * 
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Abstract. Hi tnis paper, wc> proposea numerical mettioa to approximate Trie solution or partial 

differential equations in irregular domains with no-flux boundary conditions by means of spectral 

methods. The main features of this method are its capability to deal with domains of arbitrary shape 

and its easy implementation via Fast Fourier Transform routines. We discuss several examples of 

£*\J ' practical interest and test the results against exact solutions and standard numerical methods. 
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1. Introduction. Spectral methods ^ |21 El are among the most extensively 
used methods for the discretization of spatial variables in partial differential equa- 
tions and have been shown to provide very accurate approximations of sufficiently 
smooth solutions. Because of their high-order accuracy, the use of spectral meth- 
ods has become widespread over the years in various fields, including fluid dynamics, 
quantum mechanics, heat conduction and weather prediction 01 |SJ [3 |S] . However, 
these methods have some limitations which have prevented them from being extended 
to many problems where finite-difference and finite-element methods continue to be 
used predominantly. One limitation is that the discretization of partial differential 
£> \ equations by spectral methods leads to the solution of large systems of linear or non- 

linear equations involving full matrices. Finite-difference and finite-element methods, 
on the other hand, lead to systems involving sparse matrices that can be handled by 
appropriate methods to reduce the complexity of the calculations substantially. An- 
other drawback of spectral methods is that the geometry of the problem domain must 
be simple enough to allow the use of an appropriate orthonormal basis to expand 
the full set of possible solutions to the problem. This inability to handle irregu- 
larly shaped domains is one reason why these methods have had limited use in many 
engineering problems, where finite-element methods are preferred because of their 
flexibility to describe complex geometries despite the computational costs associated 
with constructing an appropriate solution grid. Although there have been attempts 
to use spectral methods in irregular domains 9, 10 , these approaches usually involve 
k> ' either incorporating finite-element preconditioning or the use of so-called spectral el- 

ements which are similar to finite elements. We are not aware of any previous study 
where purely spectral methods, particularly those involving Fast Fourier Transforms 
(FFTs), have been used to obtain solutions in complex irregular geometries. 

In this paper we present an accurate and easy-to-use method for approximat- 
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ing the solution of partial differential equations in irregular domains with no-flux 
boundary conditions using spectral methods. The idea is based on what in dendritic 
solidification is known as the phase-field method 11 . This method is used to locate 
and track the interface between the solid and liquid states and has been applied to a 
wide variety of problems including viscous fingering |f2l I13j . crack propagation |14| 
and the tumbling of vesicles |15| . For a comprehensive review see [16| . 

In what follows we use the idea behind phase-field methods to illustrate how the 
solution of several partial differential equations can be obtained in various irregular 
and complex domains using spectral methods. Throughout the manuscript, for sim- 
plicity, we will refer to the combination of the phase-field and spectral methods as 
the spectral smoothed boundary (SSB) method. Our approach consists of two steps. 
First, the idea of the phase- field method is formalized and its convergence analyzed 
for the case of homogeneous Neumann boundary conditions. Then we discuss how 
the new formulation is useful for the direct use of spectral methods, specifically those 
based on trigonometric polynomials. This formulation makes the problem suitable 
for efficient solution using FFTs [E|. Since it is our intention that the resulting 
methodology be used in a variety of problems in engineering and applied science, we 
have concentrated on the important underlying concepts, reserving some of the more 
formal questions related to these methods for a subsequent analysis. 

2. The phase-field (smoothed boundary) method. In this work we focus 
on applying the phase-field method to partial differential equations of the form 

V{D {j) Vuj) + f(v,i, ..., ujv, t) = d t u 3 (2.1a) 

for N unknown real functions Uj defined on an irregular domain fl c R n , where 
n = 1,2,3 is the spatial dimensionality of the problem, together with appropriate 
initial conditions Uj(x,0) = Ujq(x) and subject to Neumann boundary conditions 

(H-D^Vu j ) = (2.1b) 

on dfl, where D^\x) is a family ofnxn matrices that may depend on the spatial 
variables. Equations (|2.1a|) and Ij2.1b|) include many reaction-diffusion models, such 
as those describing population dynamics or cardiac electrical activity. In here we 
will restrict the analysis to equations of the form (|2.1a() although we believe that the 
idea behind the method can be extended to many other problems involving complex 
boundaries and different types of partial differential equations. 

Instead of discretizing Eq. I|2.1a|) the smoothed boundary method relies on con- 
sidering the auxiliary problem 

V(^)r>«Vwf ) + <t>ti)f(uf\...,uf,t) = d t (^uf), (2.2) 

for the unknown functions u on an enlarged domain fl' satisfying the following 
conditions: (i) fl C fl' and (ii) dfl U dfl' = 0. The function <f>(Z' is continuous in 
fl' and has the value one inside fl and smoothly decays to zero outside fl, with £ 
identifying the width of the decay. That is, if xn is the characteristic function of the 
set fl defined as 

/ l, x e fl ,„ „n. 

*> = ! 0, xefi'-fi ^ 

then </>(*' : fl' — ► R is a regularized approximation to \n such that lim^o <j>£ = Xn ■ 
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The key idea of the Smoothed boundary method (SBM) is that when £ — > the 
solutions u of Eqs. (|2.2|) on any domain il' with arbitrary boundary conditions 



on dQ' satisfy u — ■» Uj, that means, they tend to the solution of Eqs. I|2.1a[l . 
automatically incorporating the boundary conditions (|2.1bl) . To see why, let us first 
realize that inside 17 the statement is inmediate since <fy£' — s- 1 in O as £ — » and Eq. 
<|2.2|l becomes Eq. (I2.1a|) . At the boundary we consider for simplicity the situation 
with n = 2 (the extension to n = 3 is immediate) . Assuming smoothness of dfl (which 
in this case will be a curve) and choosing any connected curve T C d£l, we define two 
families of differentiable curves T S (+) C fl'/Cl and r^-j C 12 whose ends coincide with 
those of r and which tend uniformly to T following the parameter 5. The curves T S (+) 
and r^(_) are then the boundaries of a region As whose boundary 9^4,5 = T S (+) UT$(-) 
(see Fig. l2~T|) . 

We now integrate Eq. (|2.2|) over As: 

V y,® D^Vuf ] ) + ^ f (u,t)] dx = [ f dt(^uf } )dx. (2.4) 



As- 



Using the Gauss (or Green) theorem for the first term of Eq. I|2.4|l we obtain 

I n-cf)^D i - 3) \Ju { Pdx+ f f tf>®f(u,t)dx = [ f dti^uf^dx, (2.5) 

JdA s J J As J J A s 

where § denotes a line integral over dAg. 

Now we take the limit £ — > in l|2.5|l to obtain 

lim/ n-6^D {i) Vuf l dx = -Yan j j (?) f(u,t)dx 

+ lim I [ dti^uf^dx 

= m(As) \-^f(u, t) + dti^uf)} , (2.6) 

L J J x=( 

where the last equality comes from the mean value theorem for integrals and m(As) 
is the measure of the set As. Here we assume that the solutions to Eq. 11'. 21 and its 
time derivatives are bounded so that the right-hand side of Eq. (|2.(i|) is finite. On the 
left-hand side we decompose <f a , as L + L . It is evident that lim£_,n fr n ' 

(j)^D {3) Vu { pdx = 0, since in this limit <j> ( ^ = over all T S ( +) and ft® = 1 on T S (- } . 
As we were interested in proving that the boundary conditions are satisfied, we now 
make the width of the integration region As tend to zero. Since r^(_) — ► T as 5 — > 0, 
and lim^^o in{As) = 0, we obtain 

lim lim I n ■ c# } D (j) Vu ( pdx = lim / n ■ D U) Vu ( pdx 
.5-0 «-o J aAs J s ^°Jrs^ } J 

= f n-D u} Vuf ] dx = 0. (2.7) 

Since Eq. I|2.7|) is true for any boundary segment T, we obtain the final result that in 
the limit £ — ► Eq. I|2.2|) satisfies n-D^Vuf' = for j = 1, ..., N, i.e., the boundary 
conditions. 
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Fig. 2.1. Illustration of an irregular domain Q with an example of the connected curve T and 
the domain A5 used in the proof of convergence of the method. 

From what we have shown it is clear that $7' can be any closed region containing 
Q. The idea of the smoothed boundary method is then to consider Eq. (|2.2I) for a 
small but finite £ and to discretize this problem instead of Eqs. (|2.1al) and l|2.fb|) . 
The main advantage is that one can search for the approximation m® on any enlarged 
domain O' such that D. C O' (for instance, a rectangular region). The enlarged discrete 
problem can then be solved with any proper boundary conditions on dfl', since the 
fulfilment of the boundary conditions for u on dft is guaranteed in the limit £ — > 0. In 
our case, we will use the basis of trigonometric polynomials e tkx to approximate the 
solutions; thus, we will seek an extension of the solution u of Eqs. (|2.1a() and (|2.1bjl 
that is periodic on the enlarged region Q! . 

3. The Spectral Smoothed Boundary Method. We want to discretize Eq. 
(|2.2(l on an enlarged region. As discussed earlier, we choose fi' to be a rectangular 
region containing f2 and we will expand v\^> in the basis of Cartesian products of 
trigonometric polynomials e tkxX e lk v y . 

Let us rewrite Eq. 12.21) without subscripts and superscripts as 

V(j) ■ DVu + (jN{DVu) + cf>f(u, t) = d t {(j)u). (3.1) 

Note that since <fi is located inside of the time derivative of the right term of Eq. H2.2(l , 
it is possible for the integration domain itself to evolve in time, and thus this method 
could be used to solve moving boundary problems once a coupling equation is added 
for the movement of <p. However, in this manuscript we will only deal with stationary 
integration domains; thus, dt4> = and the right side of Eq. (|3.1|l can be simplified 
as cj)d t u. Dividing Eq. (|3.1|l by </>, we get 

V log </> • DV« + V(DVu) + f(u, t) = d t u, (3.2) 

which is the equation of the smoothed boundary method that we will use to perform 
numerical simulations. 

To implement numerically any solution method for Eq. (|3.2|) . we need to make a 
specific choice for <f>^ . In practice, any method that produces a smooth characteristic 
function can be used. In the context of phase-field methods, the standard procedure 
for obtaining the values of (j>^> (which is called the "phase-field" ) is to integrate an 
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auxiliary diffusion equation of the form d t (f> — £ 2 V 2 + (20 — l)/2 — (2(j> — l) 3 /2, 
with initial conditions <f)^(t — 0) = \n, until a steady state is reached |16l I18|. 
Alternatively, since we only seek a smoothed boundary we choose to obtain <jr& from 
Xn using a convolution of the form 

(?) =Xn*G ( «\ (3.3) 

where G^ is any family of functions such that lim^o G^\x) = 6(x), where S is the 
Dirac delta function. In particular Gaussian functions of the form 

n 

G^(x) = l[exp(-xt/e) (3.4) 

fe=i 

can be chosen. In this paper, all the functions <f>^> have been obtained using this 
n-dimensional discrete convolution of xn with a Gaussian function of the form given 
by Eq. I|3.4|) . An example of the creation of <$>(& is shown in Fig. 13.11 where it can be 
seen that the width of the interface in which (j)^' changes from zero to one depends 
on the value used for £ (in fact it is of order £). 

To avoid computational difficulties for very small values of (f> we approximate 
logcj) rs log(0 + e), where e is the machine precision. Numerically, (p and (<f) + e) are 
equal up to roundoff errors, but this correction bounds the value of log(/> as <j> — > 0. 
Choosing cj) to be a periodic function to avoid Gibbs phenomenom when computing its 
derivatives forces us to choose a computational domain in which this function becomes 
small enough near <9f2'. In practice, this restriction requires us to leave a reasonable 
margin between the boundaries of the physical and the enlarged domains. We found 
that a margin of value M = 10£ is sufficient for all the simulations to be stable. 

All the spatial derivatives in Eq. (|3.2() are computed in Cartesian coordinates 
with spectral accuracy in Fourier space. If g(x) is a periodic and sufficiently smooth 
function, then its n th derivative is given by 

^ = f- 1 {{ik x ) n Hg}}> (3.5) 

where J- and T^ 1 denote the direct and inverse Fourier Transform respectively, k x are 
the wave numbers associated with each Fourier mode, and i is the imaginary unit. As 
mentioned previously, the use of this representation for u implicitly assumes periodic 
boundary conditions on <9£1'. It is significant that only Fourier Transforms are used for 
these calculations instead of differentiation matrices, thereby avoiding the generation 
and storage of these matrices and yielding more efficient codes and shorter execution 
times, especially when Fast Fourier Transforms routines are used. 

In this paper we are not concerned with designing the most efficient SSBM, but 
only to prove that such a method can be used to integrate PDEs in irregular domains. 
Thus, for time integration we use a simple second-order explicit method. In the 
particular case where all the coefficients of the diffusion tensor D are constants, we 
can write Eq. (|3.2|) in the form 

Cu + M(u,t) = d t u, (3.6) 

where Cu = V(£>Vu) is the linear term and J\f(u, t) = V log <fi ■ DVu + f(u, t) is the 
nonlinear part of the equation. Then a second-order in time operator splitting scheme 
of the form 

U(t + At) = e CAt / 2 e^ At e CAt ^ 2 U(t) (3.7) 
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Fig. 3.1. Left: Example of an irregular domain Q defined on a Cartesian grid and an enlarged 
domain Cl' . Right: Smoothing of the irregular boundary in a one- dimensional section of the domain. 
The solid line shows a small section of the characteristic function \n (with value or 1) correspond- 
ing to part of the thicker line shown in the left part of the figure. Phase-field functions <f>^> obtained 
from Xf2 for £ = 0.10, 0.05 and 0.025 are labeled by diamonds, circles and stars respectively. 



can be used to solve the equation in time |T2]- For the examples to be presented 
later, we solve the nonlinear term by a second-order (half-step) explicit method and 
integrate the linear part exactly in Fourier space by exponential differentiation, which 
reduces the stiffness of the problem considerably and allows the use of larger time 
steps. The operator splitting scheme can then be written as 



L { 



U* =T' 



,CAt/2 



Hu k }} 



+ M{U* +Af(U*,t k + At/2) ■ At/2,t k + At) ■ At 



,CAt/2 



T{U**}}. 



(3.8a) 
(3.8b) 
(3.8c) 



4. Examples of the methodology. 

4.1. The heat equation. As a first example, we will consider a simple linear 
heat equation. This first case will allow us to make a quantitative study of the errors of 
the SSB method. Specifically, we are interested in solving the following heat equation 
with sources: 



d t u = DAu-rcos(29) 



(4.1) 



in the annulus Q defined by 1 < r < 2 with homogeneous Neumann boundary con- 
ditions on dCl, d r u\ r =i = d r u\ r=2 = 0, and initial data u(r, 9, 0) = 0. The diffusion 
coefficient is taken to be constant with value D = 1. Figure |4~T1 shows one example 
of the generation of the smoothed boundary for this domain. 

Equation l|4.1|l in this geometry has an explicit steady state-solution of the form 



u s t{r,9) 



1 , 31, 8 1. 



(4.2) 
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Fig. 4.1. Left: The rectangular domain O' in which Eq. \-j.l^ is solved using the SSB method, 
with the irregular domain fi (an annulus) shown in gray. Right: Smoothed boundary function tj>^' 
given by Eq. 13.31 for £ = 0.10. 
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Fig. 4.2. Maximum absolute (left) and relative (right) errors of the numerical solution of Eq. 
\4-l I ) at time t = 6 in the annulus compared to the steady solution as a function of the parameter 
n = £/Aie. 



which can be compared with the numerical steady solution of Eq. (|4.1() (in practice, we 
stop simulations at t = 6 since by this time the numerical solutions have approximately 
reached the steady state) and obtain error estimates. Figure l4~2l shows the maximum 
absolute error E = \\u — U^\\oo and relative error e = \\u — C/^||oo/||u||oo of several 
simulations for different values of £ and grid resolutions, where the relative error is 
defined with respect to the maximum value of the analytical solution. Note that in 
general these maximum errors decrease as the thickness of the interface is reduced 
(£ — ► 0), and that in most of the simulations the relative error is less than 1%. 

Although <p&' is a continuous function, it is necessary to have a grid fine enough 
to resolve properly the boundary layers in which it quickly changes from zero to one. 
For this reason, errors are represented in Fig. 14.21 not as a function of the number of 
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Fig. 4.3. (Left) Solution of Eq. \4-l\j obtained with the SSB method at time t = 6 in the annulus 
1 < r < 2. Grid size is 540 X 540 with £ = 0.025. (Right) Spatial distribution of the absolute error 
E = \\u — U^'\\oo over the annulus for the simulation solution compared to the analytical steady 
solution. 



grid points but as a function of the parameter i] = £/A:r, which gives an idea of the 
number of points that lie in the interface. As mentioned before, it is also necessary 
to use a margin between the integration domain fi and the computational domain 17' 
for <j) to become sufficiently small on dfl'. In this particular case, 



Ax s 2{R 



N 



My 



(4.3) 



where R is the outer radius of the annulus, N is the grid resolution, and M = 10£ 
is the margin used. This expression implies that for solving Eq. I|4.1|) in the range 
7] G [1.5,5], the grid resolution varies from 90 to 300 points if £ = 0.10, from 150 to 
500 points if £ = 0.05, and equivalently from 270 to 900 points when £ = 0.025. The 
last important point concerning Fig. 14.21 is that, once the interface is properly solved 
(rj ~ 3 — 4), the error converges to an approximately constant value that depends only 
on £, so there is no reason for using excessively clustered grids. To ensure that this 
error is produced only by the spatial discretization and is not due to the order of the 
method chosen to perform the time integration, we have also run the simulations with 
a first-order explicit (Euler) time-integration method and obtained errors of the same 
order of magnitude. Figure 1431 shows the solution to Eq. I|4.1[l obtained at time t = 6 
with the SSB method. Contour lines are also included to illustrate that the no-flux 
boundary conditions at r = 1 and r = 2 are satisfied. 

The analytical solution of Eq. Q4.1[l also satisfies homogeneous Neumann bound- 
ary conditions on the quarter-annulus delimited by 1 < r < 2, < < ir/2 (see 
Fig. I4.3J1 . which allows us to use this related geometry to show how the SSB method 
performs when sharp corners are present in a given geometry. Maximum absolute and 
relative errors of the simulations for the quarter-annulus are shown in Fig. 14.41 Both 
geometries show similar good convergence properties. However, errors are slightly 
larger but of the same order of magnitude than for the full annulus due to the pres- 
ence of the sharp corners, which become slightly blunted. This can be seen in Fig. 
14.51 which shows the error distribution E = \\u — C/^- ) || 00 over the domain fi, along 
with the corresponding solution. 

In Figs. 14. 31 and 14.51 we have shown the solutions of Eq. (|4.1|l within the irregular 
geometries $7. However, the solutions U^> are calculated over the entire domain 
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Fig. 4.4. Maximum absolute (left) and relative (right) errors of the numerical solution of Eq. 
14-lf l in the quarter- annulus l<r<2, 0<$< n/2 at time t = 6 compared to the steady solution 
of the problem as a function of the parameter r\ = {/Ax. 
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Fig. 4.5. (Left) Solution to Eq. \4-l\ in the quarter annulus 1 < r < 2, < < it/2 at t = 6 
using the SSB method. Grid resolution is 500 X 500 and § = 0.025. (Right) Spatial distribution of 
absolute error E = \\u — U^Waa over the quarter-annulus for the simulation solution compared to 
the analytical steady solution. 



Q' . Figure 1-0)1 shows the solutions over Q,' in the full and quarter-annulus examples. 
While no-flux boundary conditions are implemented along <9f2, the overall solution 
has periodic boundary conditions. Note that as the solution U^' is not discontinuous 
on <9f2, our solutions never present Gibbs phenomena due to the irregular boundaries. 
An important advantage of the SSB method is that when the new formulation 
given by Eq. IJ3.2J) is used, separate equations are not written for the boundaries, as the 
solution automatically adapts to satisfy the boundary conditions on dQ, which results 
in a very simple computational implementation. Alterations to the domain geometry 
therefore can be handled straightforwardly, without generating and implementing 
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Fig. 4.6. Solution to Eq . \4-lJj ) using the SSB method over the entire domain of integration Q' 
for the annulus shown in Fig \p^ Jleft) and the quarter- annulus shown in Fig \4-5\ (right). Note that 
the periodicity along dQ' imposed by the FFT does not alter the solution in f2 or at the boundary 
d£l (shown in white) where the zero-flux boundary conditions are satisfied. 





Fig. 4.7. Solution to Eq. \4-lj at time t = 5 in a more complicated domain (left) and spatial 
distribution of the difference between the solutions obtained by the SSB method and using finite 
elements (right). Grid resolution is 400 X 400 and £ = 0.05. 



additional boundary condition equations. As an example, we present in Fig. 14.71 
the solution to Eq. i|4.1[l in a more complicated geometry that combines both polar 
and Cartesian coordinates. As we cannot obtain an analytic steady solution to the 
equation in this domain, we evolve it until time t — 5 and compare our solution to 
one obtained using the finite elements toolbox of Matlab. Since it is well known 
that spectral methods have higher accuracy than second-order finite elements, we 
obtained the finite-element solution using the finest possible grid that our machine 
could handle (31152 nodes and 61440 triangles). At this resolution, as shown in 
Fig. 14.71 the maximum difference between the solutions obtained by both methods 
is very small (s=s 10 -3 ). We expect that as the grid is refined for the finite elements 
method the solution will converge to the spectral solution, resulting in an even smaller 
difference error. 



4.2. The Allen-Cahn equation. An important partial differential equation 
which arises in the modeling of the formation and motion of phase boundaries is the 
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Allen-Cahn [20] equation: 

d t u = e 2 Au - f(u), x efl , ■. 

d„u = 0, x e an, [ ' 

where e is a small positive constant and f(u) is the derivative of a potential function 
W(u) that has two wells of equal depth. For simplicity, we will assume that W(u) = 
(u 2 — 1) 2 /4, which makes f(u) = u 3 — u. In this manner, the Allen-Cahn equation may 
be seen as a simple example of a nonlinear reaction-diffusion equation. As explained 
in Ref. 0, this equation has three fixed-point solutions, u = — 1, u = and u = 1. 
The middle state is unstable, but the states u = ±1 are attracting, and solutions tend 
to exhibit flat areas close to these values separated by interfaces that may coalesce or 
vanish on long time scales, a phenomenon known as metastability. Figure 14.81 shows 
the solution of the Allen-Cahn equation solved with Neumann boundary conditions on 
an annulus with a z-shaped hole using the SSB method, with e = 0.01. The annulus 
structure is given by 1 < r < 5 and the z-hole is formed using radii at 2, 3, and 4 
and angles in steps of 15° degrees (15, 30, 60 and 75°). For initial conditions we have 
chosen two positive and two negative Gaussian functions located in different parts of 
the domain: 

n=4 

u = ]T(-1)" +1 exp(-20((x - x t ) 2 + (y- y t ) 2 )), (4.5) 

i=i 

where 

xi = 1.5cos(7r/4), yi = 1.5sin(7r/4), 

x 2 = 4cos(tt/12), 2/2 = 4sin(7r/12), 

£3 = 4.5cos(7r/4), j/3 = 4.5sin(7r/4), 

x 4 = 4cos(llvr/24), y A = 4sin(ll7r/24). 

For comparison we also solved the equation using a second order in space finite 
difference method in polar coordinates, since despite the complexity of this geometry 
all the boundaries are parallel to the axes in polar coordinates and so implementing 
no-flux boundary conditions is straightforward. As shown in Fig. 14.81 both the SSB 
and polar finite-difference solutions are in very good agreement, except in the interface 
separating the two phases that lies just in one of the corners. These larger but still 
overall small differences between both schemes are due to the sharp transition of the 
solution between -1 and 1 and the fact that the finite difference scheme approximate 
spatial derivatives using lower-order accuracy than the spectral method. 

4.3. Reaction-diffusion equations and excitable media. Models of ex- 
citable media form another significant class of nonlinear parabolic partial differential 
equations and describe systems as diverse as chemical reactions j^ E2] , aggregation 
of amoebae in the cellular slime mold Dictyostelium |23|. calcium waves [23, and 
the electrical properties of neural |5S] and cardiac cells PHI HZ| , among others. The 
equations of excitable systems extend the Allen-Cahn equation by including one or 
more additional variables that govern growth and decay of the waves. Solutions of ex- 
citable media consist of excursions in state space from a stable rest state and a return 
to rest, with the equations describing the additional variables determining the time 
courses of excitation and recovery. In spatially extended systems, diffusive coupling 
allows excitation to propagate as nonlinear waves, and in multiple dimensions complex 
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Fig. 4.8. Solution of the Allen-Cahn equation \4-4t a * time £ = 65 with e = 10 using the SSB 
method (left) and errors when compared to the solution obtained by a second-order finite- difference 
scheme in polar coordinates (right). Grid size is 360 X 360 with £ = 0.05 for the SSB method and 
grid resolution is Ar = 0.025 and A8=0.3 degrees for the polar finite- difference solution. 



patterns can be formed, including two-dimensional spiral waves j^ E3 I2HI 1241 12~%| 
and their three-dimensional analogs, scroll waves [20 EDI- Well-known examples of 
excitable media equations include the Hodgkin-Huxley [53] model of neural cells and 
its generalized simplification, the FitzHugh-Nagumo |31j model. 

The dynamics of wave propagation in excitable media has been studied extensively 
in regular domains. However, the complex geometry inherent to some systems, such 
as the heart, often can have an essential influence on wave stability and dynamics |H2| . 
This fact, combined with the need for high-order accuracy to resolve the sharp wave 
fronts characteristic of cardiac models, should make the SSB method a useful tool for 
studying electrical waves in realistic heart geometries. Figure |4~§1 shows an example 
of a propagating wave of action potential in both an idealized (left) and a realistic 
(right) slice [33J of ventricular tissue using the SSB method and a phenomenological 
ionic cell model |3U1 132] with equations of the form 



d t u(x,t) = V ■ (DVu) - Jfi(u,v) - Jso(u) - J S i(u,w) 

d t v(x,t) = @(u c - u)(l - v)/t~[u) - Q(u - u c )v/t+ 

d t w{-x.,t) = Q(u c -u)(l -w)/t~ - <3(u-u c )w/t+ 

v 
Jfi(u,v) = 0(u - u c )(l - u)(u - u c ) 

Td 

u 1 
J so (u) = —&{u c -u)-\ 6(u - u c ) 

TO T r 



J si (u,w) 



w 
2^ 



(l+tanh[/c(u -«"')]) 



t v (u) = Q(u - u v )t v1 + 0{u v - u)t v2 



(4.6a) 
(4.6b) 
(4.6c) 

(4.6d) 
(4.6e) 
(4.6f) 
(4.6g) 



where u is the membrane potential; Jfi,J so and J S i are phenomenological currents; 
v and w are ionic gate variables; and D is the diffusion tensor (isotropic for these 
simulations, with value D — 1 cm 2 /s). In all these formulas, Q(x) is the standard 
Heaviside step function defined by Q(x) = 1 for x > and Q(x) = for x < 0, and 
the set of parameters of the model are chosen to reproduce different cellular dynamics 
measured experimentally. 
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Fig. 4.9. Propagating wave of electrical potential in an idealized (left) and a realistic (right) 
slice of ventricular tissue using the SSB method. Grid size is 400 X 400 and £ = 0.025 in both cases. 
The color code denotes tissue voltage in mV, where red corresponds to cells with higher potential 
(depolarization) and blue to cells coming back to resting state (repolarization). Following Ref. 1301 . 
parameters for these simulations in Eqs. \4-6\ are Td = 0.25, r r = 50, t s [ = 45, to = 8.3, rj" = 3.33, 
t-- = 1000, r" = 19.2, t+ = 667, r~ = 11, u c = 0.13, u v = 0.055 and uf = 0.85. 



5. Discussion and future work. In this paper we have presented a new method 
for implementing homogeneous Neumann boundary conditions using spectral meth- 
ods for several problems of general interest. The spectral smoothed boundary method 
offers several advantages over finite-difference and finite-element alternatives. Because 
ghost cells are not needed, the implementation of boundary conditions requires less 
coding than finite-difference stencils, and in addition the spatial derivatives are repre- 
sented with higher accuracy. The use of simple Cartesian grids makes the SSB method 
easier to use with multiple domain shapes than finite elements, since grid generation 
is not necessary. Furthermore, the use of FFT routines in the SSB method ensures 
efficiency and also makes extension of the method to three dimensions straightforward 
using well-established routines. Since the method is directly based on the FFT it is 
very simple to implement on high performance computers by using native parallel or 
vector FFT libraries. 

The most significant limitation of the SSB method is that the error of the method 
depends directly on the ratio of the width of the smoothed boundary £ to the spatial 
step Ax, what implies that using uniform grids the number of points in the discretiza- 
tion needs to be large. This limitation is perhaps not important for certain classes 
of problems in which the solution contains steep wavefronts or other sharp features 
that require a fine spatial resolution to correctly reproduce the dynamics of the sys- 
tem, such as electrical waves in cardiac tissue or shock waves in fluid mechanics. 
However, the adequate reduction of error in domains with irregular boundaries using 
this method may require an increase in spatial resolution of a factor of 10 or 20 in 
each direction of the mesh for problems with smooth behavior like the heat equation 
with slowly varying sources compared to what is typically needed to obtain the same 
accuracy without complex boundaries. 

Thus, an important future extension of this work is to improve the performance of 
the SSB method for problems that do not track features with sharp spatial gradients. 
For instance, if the boundary is stationary, it could be useful to use a non-uniform grid 
with extra resolution along the boundaries combined with a non-uniform fast Fourier 
transform (NFFT) to calculate the spatial derivatives. However, if the boundary 
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moves over time, it might be more efficient to use a fine spatial discretization than to 
keep track of the boundary for such problems. Other planned future work includes 
properly handling complex anisotropics in the diffusion matrices such as those found 
in cardiac muscle, examining whether the method can be used to satisfy other types of 
boundary conditions including Dirichlct and Robin, and implementing non-stationary 
boundaries. 

In conclusion, we have presented and analyzed a new numerical method which 
imposes homogeneous Neumann boundary conditions in complex geometries using 
spectral methods. We have used this method to solve different partial differential 
equations in domains with irregular boundaries and have found good agreement with 
the exact analytical solutions when such solutions can be obtained. Along with the 
overall advantage of allowing domains of different shapes to be considered with spec- 
tral methods in a very simple way, this method also offers highly accurate discretiza- 
tions of spatial derivatives, ease of implementation, straightforward extension to 3D, 
and applicability to a wide variety of equations. Moreover, SSB codes need not change 
to implement different geometries since all the information on the geometry is con- 
tained in the function <j)^\ with the additional advantage that this function is easy 
to generate and, unlike finite element methods, does not require the use of special 
software for grid generation. 
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